#include <iostream>
#include <cassert>
#include <cmath>
#include <vector>
#include <map>
#include <algorithm>
#include "euler/prime_factor.hpp"
#include "euler/prime_table.hpp"
#include "euler/gcd.hpp"
#include "euler/totient.hpp"
#include "euler.h"

BEGIN_PROBLEM(302, solve_problem_302)
	PROBLEM_TITLE("Strong Achilles numbers")
	PROBLEM_ANSWER("1170060")
	PROBLEM_DIFFICULTY(3)
	PROBLEM_FUN_LEVEL(2)
	PROBLEM_TIME_COMPLEXITY("(log(N))^(N^(1/3)*log(N))")
	PROBLEM_SPACE_COMPLEXITY("N^(1/3)")
	PROBLEM_KEYWORDS("totient")
END_PROBLEM()

#define VERBOSE 0

///////////////////////////////////////////////////////////////////////////
// Helper functions

// Displays the prime factors of n.
static void show_factors(long long n)
{
	std::cout << n << " =";
	euler::prime_factorize(n, [](long long p) { std::cout << " " << p; });
	std::cout << std::endl;
}

///////////////////////////////////////////////////////////////////////////
// Algorithm for small scale

static bool is_achilles(int n)
{
	bool powerful = true;
	int gcd_power = 0;
	euler::prime_factorize_distinct(n, [&powerful,&gcd_power](long long p, int k) {
		if (k == 1)
			powerful = false;
		gcd_power = euler::gcd(gcd_power, k);
	});
	return powerful && (gcd_power==1);
}

static bool is_strong_achilles(int n)
{
	return is_achilles(n) && is_achilles(euler::totient(n));
}

static int verify(int N)
{
	int count = 0;
	for (int n = 2; n <= N; n++)
	{
		if (is_strong_achilles(n))
		{
#if VERBOSE
			show_factors(n);
#endif
			++count;
		}
	}
	return count;
}

///////////////////////////////////////////////////////////////////////////
// Algorithm for large scale

// For each prime factor p, do a prime-factorization of (p-1).
// The list of all factors (p-1) are stored consecutively in _factors_.
// The start position of the i-th factorization is stored in _factor_pos_.
static void factorize_p_minus_one(
	const std::vector<int> &primes,
	std::vector<std::pair<int,int>> &factors,
	std::vector<int> &factor_pos)
{
	factors.resize(0);
	factor_pos.resize(primes.size()+1);

	for (int i = 1; i < (int)primes.size(); i++)
	{
		int t = primes[i] - 1;
		for (int j = 0; ; j++)
		{
			int pp = primes[j];
			if (pp*pp > t)
				break;
			int kk = 0;
			while (t % pp == 0)
			{
				t /= primes[j];
				kk++;
			}
			if (kk > 0)
				factors.push_back(std::pair<int,int>(j, kk));
		}
		if (t > 1)
		{
			int j = (int)(std::lower_bound(primes.cbegin(), primes.cend(), t) - primes.cbegin());
			factors.push_back(std::pair<int,int>(j, 1));
		}
		factor_pos[i+1] = (int)factors.size();
	}
}

static unsigned long long square(int x)
{
	return (unsigned long long)x*x;
}

static unsigned long long cube(int x)
{
	return (unsigned long long)x*x*x;
}

static int square_root(long long n)
{
	int r = (int)sqrt((long double)n)+1;
	if ((long long)r*r > n)
		r--;
	return r;
}

static int cube_root(long long n)
{
	int r = (int)pow((long double)n, (long double)1.0/3.0)+1;
	if ((long long)r*r*r > n)
		r--;	
#if 0
	if (!((long long)r*r*r <= n && (long long)(r+1)*(r+1)*(r+1) > n))
		abort();
#endif
	return r;
}

class problem_302
{
private:
	typedef unsigned long long uint64_t;

	long long N;				// max n
	std::vector<int> P;			// prime factor base
	std::vector<uint64_t> P2;	// p^2
	std::vector<uint64_t> P3;	// p^3
	std::vector<std::pair<int,int>> factors; // factorization of (p-1)
	std::vector<int> factor_pos;
	unsigned char GCD[64][64];  // gcd table

	int search(uint64_t B, uint64_t n, int i, int gcd_K, int gcd_L, std::map<int,int> &L) const;

public:
	problem_302(long long max_n);
	long long solve() const;
};

problem_302::problem_302(long long max_n)
{
	N = max_n;
	euler::prime_table<int> ptable(cube_root(N/4));
	P.assign(ptable.begin(), ptable.end());
	factorize_p_minus_one(P, factors, factor_pos);

	// Cache p^3.
	P3.resize(P.size());
	for (int i = 0; i < (int)P.size(); i++)
	{
		P3[i] = cube(P[i]);
	}

	// Cache GCD results.
	for (int a = 0; a < 64; a++)
	{
		for (int b = 0; b < 64; b++)
		{
			GCD[a][b] = euler::gcd(a, b);
		}
	}
}

int problem_302::search(
	uint64_t bound,
	uint64_t n,
	int i,
	int gcd_K, int gcd_L,
	std::map<int,int> &L) const
{
	assert(i > 0);

	int count = 0;

	// Find the largest p_i that can possibly become a factor of n.
	int imax = (int)(std::lower_bound(P3.cbegin(), P3.cbegin()+i, bound+1) - P3.cbegin() - 1);
	int p2 = square_root(bound);
	for (auto it = L.crbegin(); it != L.crend(); ++it)
	{
		int ii = it->first;
		if (ii >= i)
			continue;
		if (ii <= imax)
			break;

		int pp = P[ii];
		if (pp <= p2)
		{
			imax = ii;
			break;
		}

		int kk = it->second;
		if (kk == 1) // phi(n) contains the factor with power 1, but n can't
			return 0;
		gcd_L = GCD[gcd_L][kk];
	}

	// Assume each p_i to be the largest factor in n.
	for (i = imax; i >= 0; i--)
	{
		// p[i] can be included if p_i^k_i <= Bound.
		int k_i;
		uint64_t prod;
		auto it_L = L.find(i);
		int L_i = (it_L != L.end())? it_L->second : 0;
		if (L_i == 0) // this can be optimize by adding a ptr for L
		{
			k_i = 3;
			prod = cube(P[i]);
		}
		else
		{
			k_i = 2;
			prod = square(P[i]);
		}
		if (prod <= bound)
		{
			uint64_t B = bound / prod;

			// Now we want to include the prime p[i] in n. Then 
			// we need to factorize (p[i]-1) and update the vector L.
			// During this process, if we find the largest L[i]=1,
			// then we know that we must have a prime factor q>=p[i]
			// where q^2 exists.
			int p1 = 1;
			for (int j = factor_pos[i]; j < factor_pos[i+1]; j++)
			{
				 if ((L[factors[j].first] += factors[j].second) == 1)
					 p1 = P[factors[j].first];
			}
			uint64_t LB = square(p1);

			// Try each admissible power of p[i].
			for (int l_i = L_i + k_i - 1; B >= LB; B /= P[i], prod *= P[i], k_i++, l_i++)
			{
				int _gcd_k = GCD[gcd_K][k_i];
				int _gcd_l = GCD[gcd_L][l_i];
				if (i == 0)
				{
					if (_gcd_k == 1 && _gcd_l == 1)
					{
#if VERBOSE
						std::cout << "Found: " << (n*prod) << std::endl;
#endif
						count++;
					}
				}
				else
				{
					count += search(B, n*prod, i, _gcd_k, _gcd_l, L);
				}
			}

			// Restore the values of L[i].
			for (int j = factor_pos[i]; j < factor_pos[i+1]; j++)
			{
				if ((L[factors[j].first] -= factors[j].second) == 0)
					L.erase(factors[j].first);
			}
		}
		
		// p[i] has to be included if L[i] = 1.
		if (L_i == 1)
			return count;

		// Update gcd_L.
		gcd_L = GCD[gcd_L][L_i];
	}

	// If no more prime factors could be possibly included in n,
	// we return from here.
#if VERBOSE
	if (gcd_K == 1 && gcd_L == 1)
	{
		std::cout << "Found: " << n << std::endl;
	}
#endif
	count += (gcd_K == 1 && gcd_L == 1)? 1 : 0;
	return count;
}

long long problem_302::solve() const
{
	std::map<int,int> L;
	return search(N, 1, (int)P.size(), 0, 0, L);
}

// Three optimizations:
// 1. remove the recursion of the last level
// 2. use a map to store the power presentation; then find the max prime as 
//    the max of:
//    - p_i where L[i]=0 and p_i^3<bound
//    - p_i where L[i]>0 and p_i^2<bound
//    The p_i found this way must be larger than p_j where L[j]=1
// 3. use static allocation for the map
static void solve_problem_302()
{
	//const long long N = 10000;
	//const long long N = 100000000;
	const long long N = 1000000000000000000LL; // 10^18

	//std::cout << "Verify count: " << verify(N) << std::endl;
	std::cout << problem_302(N).solve() << std::endl;
}
